package com.asa.numanaly;

import com.asa.HanShu;
import com.asa.utils.Calculate;

/**
 * 逼近论 p431 函数比较拟合，这个和之前C的区别是，C是一定会经过所有的点，用函数来描述图像
 * 
 * @author Administrator
 *
 */
public class H {

	/**
	 * Pade有理逼近
	 * 花了我很多时间
	 * 把q0、p0都输出了，但实际使用不需要这么多数据
	 * 
	 * (p[0]+p[1]*x+p[2]*x^2+p[3]*x^3)/(q[0]+q[1]*x+q[2]*x^2)
	 * 
	 * @return
	 */
	public static double[][] bijing1(int m, int n, double[] y) {
		int N = m + n;
		double[] a = new double[N + 1];
		double[] q = new double[m + 1];
		double[] p = new double[n + 1];
		double[][] b = new double[N + 1][N + 2];

		for (int i = 0; i <= N; i++) {
			a[i] = Math.pow(-1,i)*Math.pow(y[i],i) / Calculate.jiecheng(i);
		}
		System.out.println("a的值");
		for (int i = 0; i < a.length; i++) {
			System.out.println(a[i]);
		}
		System.out.println("====================");
		
		

		q[0] = 1;
		p[0] = 1;

		for (int i = 0; i <= N; i++) {

			for (int j = 1; j <= i - 1; j++) {
				if (j <= n) {
					b[i][j] = 0;

				}
			}

			if (i <= n) {
				b[i][i] = 1;
			}

			for (int j = i + 1; j <= N; j++) {
				b[i][j] = 0;
			}
			for (int j = 1; j <= i; j++) {
				if (j <= m) {
					b[i][n + j] = -a[i - j];
				}
			}

			for (int j = n + i + 1; j <= N; j++) {
				b[i][j] = 0;
			}

			b[i][N + 1] = a[i];

		}
		System.out.println("b的初始化");
		for (int i = 0; i < b.length; i++) {
			for (int j = 0; j < b[i].length; j++) {
				System.out.print(b[i][j]+"      ");
			}
			System.out.println();
		}
		
		
		double[][] fca = new double[N][N];
		double[] fcb = new double[N];

		for (int i = 0; i < fca.length; i++) {
			for (int j = 0; j < fca.length; j++) {
				
				fca[i][j] = b[i+1][j+1];
				
				
			}
		}
		for (int i = 0; i < fcb.length; i++) {
			fcb[i] = a[i+1];
		}
		System.out.println("00000000000000000000000");

		
		double[] fangchengzu1 = F.fangchengzu1(fca, fcb);
		
		for (int i = 0; i < fangchengzu1.length; i++) {
			System.out.println(fangchengzu1[i]);
		}
		System.out.println("00000000000000000000000");

		for (int i = 0; i < n; i++) {
			System.out.println(fangchengzu1[i]);
			p[i+1] = fangchengzu1[i];

		}
		System.out.println("00000000000000000000000");

		for (int i = n; i < N; i++) {
			System.out.println(fangchengzu1[i]);
			q[i-n] = fangchengzu1[i];
		}
		
		double[][] result = new double[2][];
		result[0] = p;
		result[1] = q;
		return result;
		

	}
	
	
	/**
	 * Chebushev有理逼近
	 * 格式化成矩阵没有大问题，但矩阵求解不对
	 * 有偏差，但还算是同位上的偏差，应该是由于积分导致的误差
	 * 这里面还有Ti(x)的概念，要使用一定要看一下文档，这个算法论文中提示有比其跟精准的	Remes算法
	 * @param m
	 * @param n
	 * @param y
	 * @return
	 */
	public static double[][] bijing2(int m, int n, double[] y) {
		int N = m + n;
		double[] a = new double[y.length];
		double[] q = new double[m + 1];
		double[] p = new double[n + 2];
		double[][] b = new double[N + 1][N + 2];

		for (int i = 0; i < a.length; i++) {
			a[i] = y[i];
		}
		System.out.println("a的值");
		for (int i = 0; i < a.length; i++) {
			System.out.println(a[i]);
		}
		System.out.println("====================");
		
		

		q[0] = 1;

		for (int i = 0; i <= N; i++) {

			for (int j = 0; j <= i; j++) {
				if (j <= n) {
					b[i][j] = 0;

				}
			}

			if (i <= n) {
				b[i][i] = 1;
			}

			for (int j = i + 1; j <= n; j++) {
				b[i][j] = 0;
			}
			for (int j = n+1; j <= N; j++) {
				if (i!=0) {
					b[i][j] = -(a[i+j-n]+a[Math.abs(i-j+n)])/2;
				}else {
					b[i][j] = -a[j-n]/2;

				}
			}
			
			
			if (i!=0) {
				b[i][N+1] = a[i];
				
			}else {
				b[i][N+1] = a[i]/2;

			}
			


		}
		System.out.println("b的初始化");
		for (int i = 0; i < b.length; i++) {
			for (int j = 0; j < b[i].length; j++) {
				System.out.print(b[i][j]+"      ");
			}
			System.out.println();
		}
		
		
		double[][] fca = new double[N+1][N+1];
		double[] fcb = new double[N+1];
		System.out.println("初始化A的部分：");
		for (int i = 0; i < fca.length; i++) {
			for (int j = 0; j < fca.length; j++) {
				
				fca[i][j] = b[i][j];
				System.out.print(b[i][j]+"      ");

				
			}
			System.out.println();
		}
		System.out.println("初始化B的部分：");

		for (int i = 0; i < fcb.length; i++) {
			fcb[i] = b[i][N+1];
//			if (i==fcb.length-1) {
//				fcb[i] = 0.000543;
//			}
			System.out.println(b[i][N+1]);
		}
		System.out.println("00000000000000000000000");

		
		double[] fangchengzu1 = F.fangchengzu1(fca, fcb);
		
		for (int i = 0; i < fangchengzu1.length; i++) {
			System.out.println(fangchengzu1[i]);
		}
		System.out.println("ppppppppppppppppppppp");

		for (int i = 0; i <= n; i++) {
			System.out.println(fangchengzu1[i]);
			p[i+1] = fangchengzu1[i];

		}
		System.out.println("qqqqqqqqqqqqqqqqqqqqqqqqqqq");

		for (int i = n+1; i < fangchengzu1.length; i++) {
			System.out.println(fangchengzu1[i]);
			q[i-n] = fangchengzu1[i];
		}
		
		double[][] result = new double[2][];
		result[0] = p;
		result[1] = q;
		return result;
		

	}
	
	
	
	
	/**
	 * 傅里叶变换
	 * 用三角函数来模拟函数
	 * p477		书上面没有算法，但有计算公式
	 * @param m	传入是5 那么a计算出来有4 b计算出来有2		在我这里a的首尾提出来，中间1、2和b配对
	 * @param p
	 * @param y
	 * @return 第一行是a第二行是b
	 * 
	 * s(z) = a[0]/2+a[m-2]*cos((m-2)*z)+连加符合（k=1到m-3）(a[k]*cos(k*z)+b[k]*sin(k*z))
	 * 其中z=PI*(x-1)
	 */
	public static double[][] bijing3_0(int m,double[] x, double[] y) {
		
		
		double[] z = new double[x.length];
		
		for (int i = 0; i < z.length; i++) {
			z[i] = Math.PI*(x[i]-1);
		}
		
		double[] a = new double[m-1];
		double[] b = new double[m-2];
		
		
		for (int k = 0; k < a.length; k++) {
			
			double asa = 0;
			for (int j = 0; j < x.length; j++) {
				asa += y[j]*Math.cos(k*z[j]);
			}
			
			a[k] = asa/5;
			System.out.println(a[k]);
		}
		
		System.out.println("==========================");
		for (int k = 1; k < b.length; k++) {
			double asa = 0;
			for (int j = 0; j < x.length; j++) {
				asa += y[j]*Math.sin(k*z[j]);
			}
			b[k] = asa/5;
			System.out.println(b[k]);
		}
		
		
		
		
		
		return null;
		

	
	
	}
	
	
	
	
	
	/**
	 * 快速Fourier变换
	 * 
	 * 快速傅里叶变换
	 * 算法中有不清晰的地方，我要再找算法才行
	 * 这个的原理是把三角函数通过欧拉方程，放到了复数域，重点是复数域
	 * 参考文献【Ham】
	 * 参考文献【AHU,pp.252-269】
	 * @param m
	 * @param n
	 * @param y
	 * @return
	 */
//	public static double[][] bijing3(int m, double p, double[] y) {
//		
//		
//		int M = m;
//		double q = p;
//		double[] asa = new double[2*m];
//		double[] c = new double[2*m];
//
//		for (int i = 0; i < 2*m; i++) {
//			asa[i] = Math.pow(Math.E, Math.PI*i/m);
//		}
//		
//		for (int j = 0; j < 2*m; j++) {
//			c[j]=y[j];
//		}
//		
//		double[] asb = new double[M+1];
//		for (int j = 1; j <= M; j++) {
//			asb[j] = asa[j];
//			asb[j+M] = -asb[j];
//		}
//		int K = 0;
//		asb[0]=1;
//		for (int L = 1; L <= p+1; L++) {
//			
//			while (K<2*m-1) {
//				for (int j = 1; j <= M; j++) {
//					 
//					K = k[p]
//					
//					
//					
//				}
//				
//			}
//			
//			
//		}
//		
//		
//		
//		
//		
//		return null;
//		
//		
//		
//	}
//	
	
	
	
	
	
	
	
	
	
	
	/**
	 * 一般性方法,类似于待定系数，函数拟合、
	 * 
	 * 这种形式
	 * y = −1.04103 ln x − 1.26132 cos x + 0.03073ex
	 * 
	 * @param args
	 */
	public static double[] Norm(double[] x,double[] y) {
		
		
		double[][] a = new double[3][3];
		double[] b = new double[3];
		
		
		a[0][0] = oneone(x);
		a[0][1] = onetwo(x);
		a[0][2] = onethree(x);
		a[1][0] = twoone(x);
		a[1][1] = twotwo(x);
		a[1][2] = twothree(x);
		a[2][0] = threeone(x);
		a[2][1] = threetwo(x);
		a[2][2] = threethree(x);
		
		
		b[0] = one(x, y);
		b[1] = two(x, y);
		b[2] = three(x, y);
		
		
		for (int i = 0; i < a.length; i++) {
			for (int j = 0; j < a.length; j++) {
				System.out.print(a[i][j]+"  ");
			}
			System.out.println();
		}
		for (int i = 0; i < b.length; i++) {
			System.out.print(b[i]+"  ");
		}
		System.out.println("解这个方程组再回调将参数回填就可得函数拟合");
		
		
		double[] fangchengzu1 = F.fangchengzu1(a, b);

		
		for (int i = 0; i < fangchengzu1.length; i++) {
			System.out.println(fangchengzu1[i]);
		}
		
		
		return fangchengzu1;
		
		
	}
	
	public static HanShu Norm2hanshu(final double[] cangshu) {
		
		if (cangshu==null||cangshu.length!=3) {
			throw new RuntimeException("传入参数不能为null或者长度不能错误");
		}

		HanShu hanShu = new HanShu() {
			
			@Override
			public double hanshu(double t, double x) {
				// TODO Auto-generated method stub
				return 0;
			}
			
			@Override
			public double hanshu(double x) {
				// TODO y = −1.04103 ln x − 1.26132 cos x + 0.03073ex
				return cangshu[0]*Math.log(x) + cangshu[1]*Math.cos(x) + cangshu[2]*Math.exp(x) ;
			}
		};
		
		
		return hanShu;
	}

	
	
	
	
	
	
	
	
	
	
	public static void main(String[] args) {
//		test1();
		double[] x = { 0.24, 0.65, 0.95, 1.24, 1.73, 2.01, 2.23, 2.52, 2.77,
				2.99 };
		double[] y = { 0.23, -0.26, -1.10, -0.45, 0.27, 0.10, -0.29, 0.24,
				0.56, 1.00 };
		Norm(x,y);
		
		
		
		
		
		
		
	}


	private static void test1() {
		HanShu hanShu = new HanShu() {
			
			@Override
			public double hanshu(double t, double x) {
				// TODO Auto-generated method stub
				return 0;
			}
			
			@Override
			public double hanshu(double x) {
				// TODO Auto-generated method stub
				return Math.pow(x, 4)-3*Math.pow(x, 3)+2*Math.pow(x, 2)-Math.tan(x*(x-2));
			}
		};
		
//		double[] x = {0.125,0.375,0.625,0.875,1.125,1.375,1.625,1.875};
//		double[] y = {0.26440,0.84081,1.36450,1.61282,1.36672,0.71697,0.17909,-0.14576};
		
		int n = 10;
		double[] x = new double[n];
		double[] y = new double[n];
		for (int i = 0; i < n; i++) {
			
			x[i] = (double)i/5.0;
			y[i] = hanShu.hanshu(x[i]);
		}
		
		
		bijing3_0(5, x, y);
	}
	
	
	
	
	

//	public static void main(String[] args) {
//
//		final HanShu hanShu = new HanShu() {
//			
//			@Override
//			public double hanshu(double t, double x) {
//				// TODO Auto-generated method stub
//				return 0;
//			}
//			
//			@Override
//			public double hanshu(double x) {
//				// TODO Auto-generated method stub
//				
//				return Math.pow(Math.E, -x);
//			}
//		};
//		
//		int m = 2;
//		int n = 3;
//		double[] x = new double[6];
//		double[] y = new double[6];
//		double b = 0.2;
//		for (int i = 1; i < y.length + 1; i++) {
//			x[i - 1] = b * i;
////			y[i - 1] = Math.pow(Math.E, -x[i - 1]);
//			y[i - 1] = hanShu.hanshu(0);
//			System.out.println(x[i - 1]+"           "+y[i - 1]);
//		}
//
////		bijing1(m, n, y);
//		
//		
//		
//		
//		HanShu hanShuy = new HanShu() {
//			
//			@Override
//			public double hanshu(double k, double x) {
//				// TODO Auto-generated method stub
//				return hanShu.hanshu(Math.cos(x))*Math.cos(k*x);
//			}
//			
//			@Override
//			public double hanshu(double x) {
//				// TODO Auto-generated method stub
//				
//				return hanShu.hanshu(Math.cos(x));//*Math.cos(x);
//			}
//		};
//		
//		
//		
//		double[] fy = new double[2*n+m+1];
//	
//		fy[0] = Calculate.hanshujifen2(0.0, Math.PI, hanShuy)*2/Math.PI;//可以考虑用D里面的积分求解，但那个没有测试过，这个粗糙点
//				
//		for (int k = 1; k < fy.length; k++) {
//			HanShu kfunction = kfunction(hanShuy, k);
//			fy[k] = Calculate.hanshujifen(0.0, Math.PI, kfunction)*2/Math.PI;;
//			
//		}
//		System.out.println("输出fy，这些奇奇怪怪的导数");
//		for (int i = 0; i < fy.length; i++) {
//			System.out.println(fy[i]);
//		}
//		
//		
//		bijing2(m, n, fy);
//
//	}
//
//
//	private static HanShu kfunction(final HanShu hanShu,final double k) {
//		System.out.println(k+"      kkkkkkkkkkkkkkkkkkkk");
//
//		HanShu hanShuyk = new HanShu() {
//			@Override
//			public double hanshu(double t, double x) {
//				// TODO Auto-generated method stub
//				return 0;
//			}
//			
//			@Override
//			public double hanshu(double x) {
//				// TODO Auto-generated method stub
//				return hanShu.hanshu(x)*Math.cos(k*x);
//			}
//
//
//			
//		};
//		return hanShuyk;
//	}
	
	
	
	
	
	
	
	
	
	
	public static double oneone(double[] x) {
		
		double all = 0;
		for (int i = 0; i < x.length; i++) {
			all = all+Math.log(x[i])*Math.log(x[i]);
		}
		
		
		return all;
	}

	public static double onetwo(double[] x) {
		
		double all = 0;
		for (int i = 0; i < x.length; i++) {
			all = all+Math.log(x[i])*Math.cos(x[i]);
		}
		
		return all;
	}

	public static double onethree(double[] x) {
		double all = 0;
		for (int i = 0; i < x.length; i++) {
			all = all+Math.log(x[i])*Math.exp(x[i]);
		}
		
		return all;
	}

	
	public static double twoone(double[] x) {

		double all = 0;
		for (int i = 0; i < x.length; i++) {
			all = all+Math.log(x[i])*Math.cos(x[i]);
		}
		
		return all;
	}

	public static double twotwo(double[] x) {
		double all = 0;
		for (int i = 0; i < x.length; i++) {
			all = all+Math.cos(x[i])*Math.cos(x[i]);
		}
		
		return all;
	}
	public static double twothree(double[] x) {
		double all = 0;
		for (int i = 0; i < x.length; i++) {
			all = all+Math.cos(x[i])*Math.exp(x[i]);
		}
		
		return all;
	}


	public static double threeone(double[] x) {
		double all = 0;
		for (int i = 0; i < x.length; i++) {
			all = all+Math.log(x[i])*Math.exp(x[i]);
		}
		
		return all;
	}

	public static double threetwo(double[] x) {
		double all = 0;
		for (int i = 0; i < x.length; i++) {
			all = all+Math.cos(x[i])*Math.exp(x[i]);
		}
		
		return all;
	}

	public static double threethree(double[] x) {
		double all = 0;
		for (int i = 0; i < x.length; i++) {
			all = all+Math.exp(x[i])*Math.exp(x[i]);
		}
		
		return all;
	}
	
	public static double one(double[] x,double[] y) {
		double all = 0;
		for (int i = 0; i < x.length; i++) {
			all = all+y[i]*Math.log(x[i]);
		}
		
		return all;
	}
	public static double two(double[] x,double[] y) {
		double all = 0;
		for (int i = 0; i < x.length; i++) {
			all = all+y[i]*Math.cos(x[i]);
		}
		
		return all;
	}
	public static double three(double[] x,double[] y) {
		double all = 0;
		for (int i = 0; i < x.length; i++) {
			all = all+y[i]*Math.exp(x[i]);
		}
		
		return all;
	}
	
	
	
	

}
